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1.  Introduction 


The  characterization  of  the  harmonic  scattering  response  of  circuit-component-loaded  scatterers 
is  an  important  problem  in  both  sensing  and  electromagnetic  interference  applications.  A  circuit- 
component-loaded  scatterer  is  defined  herein  as  a  scatterer  that  contains  circuit  lumped-elements 
(e.g.,  resistors,  inductors,  capacitors)  and/or  radio  frequency  (RF)  components  (e.g.,  transmission 
lines,  antennas,  distributed  filters,  amplifiers)  within  its  structure.  In  contrast  to  the  standard 
radar  scattering  problem,  which  necessitates  the  analysis  of  the  radio  responses  of  objects  solely 
based  on  their  dielectric  compositions  and  physical  geometries — and  which  has  been  thoroughly 
treated  by  numerous  studies,  the  generalized  scattering  characterization  of  circuit-component- 
loaded  objects  has  thus  far  garnered  relatively  less  attention  in  open  literature — the  interesting 
applications  to  which  this  research  area  can  be  applied  notwithstanding.  From  the  sensing 
perspective,  a  circuit-component-loaded  scatterer — when  properly  excited — could  potentially  re¬ 
radiate  unique  features,  enabling  it  to  be  differentiated  from  an  ordinary  scatterer:  for  example, 
the  antenna  and  filter  at  the  front-end  of  an  RF  device  can  modulate  the  radar  backscattering 
response,  facilitating  opportunities  for  classification  and  detection;  the  nonlinear  components 
within  the  scatterer  can  also  generate  harmonic  signals  and  re-emission,  which  can  be  exploited 
to  extract  the  object  embedded  in  linear  clutter.  As  related  to  the  design  of  communication 
systems,  circuit-element  modulated  scattering  plays  a  role  in  understanding  how  inter-coupling 
effects  among  different  radios  influence  transceiver  performance  and  signal  fidelity.  In  all  the 
aforementioned  applications,  an  end-to-end  (or  transmitter  terminal-to-receiver  terminal), 
system-level  simulation  capability  is  needed  to  attain  a  fundamental  understanding  of  the  signal 
interactions.  Accordingly,  in  this  study,  a  systematic  strategy  for  analyzing  the  linear  and 
nonlinear  responses  of  the  circuit-component-loaded  scatterer  is  developed,  with  specific 
relevance  to  active  RF  sensing,  detection,  and  localization. 

The  formidable  scope  of  the  current  problem  inherently  spans  both  the  electromagnetic  and 
circuit  domains,  and  a  full-wave  analysis  strategy  must  track  the  signals  or  waves  across  a 
multiscale  simulation  domain  consisting  of  effects  occurring  at  the  propagation,  RF  component, 
and  circuit  levels.  As  shown  in  figures  1  and  2,  in  the  external  problem,  the  propagation  and 
scattering  interactions  of  the  radiation  (or  current)  sources  with  the  environment  determine  the 
incident  fields  at  the  scatterer  as  well  as  the  re-radiated  fields  subsequently  received  at  the 
receiver.  In  the  internal  problem,  the  coupling  of  energy  into  the  scatterer  and  the  generation  of 
harmonics  are  complicated  functions  of  the  circuit  elements  and  RF  components.  Full-wave 
electromagnetic  solvers — for  instance,  the  method-of-moments  (MoM)  and  finite-difference 
time-domain  (FDTD) — can  be  used  to  characterize  RF  components  such  as  antennas  and 
distributed  filters.  The  treatment  of  active,  passive,  linear,  and  nonlinear  lumped-elements  also 
has  been  considered  in  the  FDTD  framework  (1-4).  In  reference  5,  a  harmonic  balance 
procedure — in  conjunction  with  MoM — is  proposed  to  analyze  the  radar  cross  section  of  simple 
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wire  scatterers  with  nonlinear  loads.  To  predict  the  microwave  nonlinear  behaviors  of  high 
temperature  superconducting  filters,  an  iterative  MoM  approach  is  developed  in  reference  6:  the 
linear  response  is  calculated  in  the  frequency  domain  using  the  electric  field  integral  equation 
and  the  nonlinear  one  is  solved  in  the  time  domain.  Works  of  direct  relevance  to  RF  transceiver 
design  also  have  been  undertaken  by  various  researchers:  the  scattering  of  antennas  connected  to 
nonlinear  loads  is  addressed  in  7-77;  harmonic  generation  by  reactively-tunable  antennas  is 
investigated  via  a  full- wave  approach  in  1 8  as  well  as  a  mixture  of  full- wave  and  equivalent- 
circuit  methods  in  19.  As  related  to  electromagnetic  inference  and  compatibility  studies,  the 
prediction  of  emissions  from  printed  circuit  boards  (PCBs)  using  numerical  methods  (e.g.,  partial 
element  equivalent  circuit  method,  finite  integration  technique,  FDTD)  has  been  discussed  in 
works  by,  for  example,  20-23 — albeit  the  focus  of  these  studies  is  not  on  nonlinear  effects.  As 
alluded  to  above,  although  there  exists  a  panoply  of  works — in  both  the  electromagnetic  and 
circuit  domains — that  are  potentially  applicable,  either  directly  or  indirectly,  to  the  treatment  of 
loaded  scatterers,  a  complete  system-level  analysis  of  the  harmonic  scattering  problem — that  is, 
the  tracing  of  the  signal  interaction  from  the  transmitter  terminal,  across  an  environment,  into  an 
RF  circuit  network,  and  back  to  the  receiver — has  not  been  considered  in  detail. 


Figure  1 .  Standoff  sensing  of  circuit-component-loaded  scatterer  buried  in  the  ground. 


2 


Figure  2.  Standoff,  through- wall  sensing  of  circuit-component-loaded  scatterer. 

(A  flat  roof  structure  is  assumed  but  not  shown.) 

In  view  of  the  current  lack  of  a  comprehensive,  generalized  solver  for  studying  the  stated 
multiscale  radar  scattering  effects,  a  hybrid  simulation  approach  is  presented  in  this  work — with 
an  emphasis  on  the  characterization  of  the  scattering  properties  of  circuit-component-loaded 
scatterers  in  complex  environments.  The  framework  of  the  analysis  involves  a  combination  of 
electromagnetic  and  circuit  simulations:  the  linear  (or  ordinary)  response  is  determined  in  the 
electromagnetic  domain,  whereas  the  nonlinear  response  is  deduced  by  computing  the  currents 
and  voltages  of  an  equivalent  circuit  network  in  a  Simulation  Program  with  Integrated  Circuit 
Emphasis  (SPICE)  solver,  and  the  interaction  of  the  electromagnetic  fields  at  the  propagation 
level  is  characterized  by  asymptotic  techniques.  The  study  is  organized  as  follows:  in  section  2, 
the  development  of  the  individual  components  of  the  simulation  framework  is  detailed;  in  section 
3,  the  application  of  the  solver  is  illustrated  for  two  practical  scenarios  of  interest — viz.,  in  the 
first  case,  the  scatterer  is  in  the  presence  of  a  dielectric  ground;  in  the  second  case,  the  scatterer  is 
in  a  room  environment;  and  in  section  4,  a  summary  of  the  work  is  presented. 


2.  Formulation  of  Hybrid  Simulation  Framework 


The  linear  and  nonlinear  responses  of  a  circuit-component-loaded  structure  (excited  by  a  time- 
harmonic  field)  are  calculated  with  a  solution  approach  encompassing  both  the  electromagnetic 
and  circuit  domains.  As  depicted  in  figure  3,  the  solver  is  separable  into  the  following  five  steps: 
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1.  The  5-parameters  and  the  current  densities  of  the  individual  RF  components  (antennas, 
distributed  filters,  PCB  structure,  etc.)  are  characterized  using  a  full-wave  approach — in 
this  case,  a  MoM  solver  is  employed.  The  initial  incident  field  at  the  scatterer  is  deduced  by 
applying  an  asymptotic  field  propagator  to  the  radiation  source  (or  transmitting  antenna). 
Note  that  to  accomplish  this  step,  the  nonlinear  circuit  elements  first  have  to  be  extracted 
from  the  circuit,  and — for  the  calculation  of  the  current  densities — the  RF  components 
must  be  terminated  with  appropriate  loads  as  determined  by  the  actual  circuit  layout. 

2.  Using  the  electromagnetic  simulation  results,  wideband  equivalent  circuit  sub-networks  are 
synthesized  for  the  RF  components.  These  equivalent  sub-networks  are  combined  with  the 
extracted  nonlinear  components  (and  other  lumped  and  discrete  elements)  from  the  original 
circuit  layout  to  form  an  overall  equivalent  network  representation  for  the  RF  structure. 

3.  The  stimulus  (voltage  or  current  source)  for  the  derived  equivalent  network  is  calculated 
based  on  the  linear  electromagnetic  response  of  the  scatterer. 

4.  The  equivalent  network  is  placed  in  a  SPICE  solver,  and  simulation  is  performed  to  obtain 
the  effective  radiated  power  at  all  the  harmonics  of  interest. 

5.  Finally,  the  re-radiation  for  each  spectral  component  is  propagated  to  the  receiver  location 
using  asymptotic  routines.  It  should  be  mentioned  that  the  first-order  solution  (for  the 
fundamental  frequency  scattering  response)  can  be  approximated  by  applying  steps  1  and  5 
only — specifically,  a  linear  MoM  approach  can  be  used  to  simulate  the  entire  RF  structure, 
with  the  discrete  circuit  elements  appearing  only  as  equivalent  impedance  loads  at  the 
structure  ports. 

Electromagnetic  Circuit 

Establish  current  density  on  structure 
using  linear-MoM  (with  lumped 
elements/equivalent  load  impedances) 


Synthesize  equivalent  network  for  original 
circuit  layout  (lumped  elements  for  discrete 
components,  vector-fitting  with  S-parameters 
for  RF  components) 


Deduce  stimulus  for 
equivalent  network 

l 

Simulate  equivalent  network  using 
circuit  solver  (PSPICE,  ADS) 


Calculate  radiated/backscattered 
field  using  MoM/asymptotics 


Figure  3.  Framework  of  hybrid  solver. 
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The  development  of  the  major  components  of  the  simulation  methodology  is  discussed  below. 

2.1  Electromagnetic  Solver 

The  linear  MoM  algorithm  for  calculating  the  RF  response  of  a  structure  is  based  on  an  approach 
originally  established  by  Liao  (24)  for  analyzing  grounded  antenna  performance.  The  details  of 
the  electromagnetic  solver  are  described  below;  for  further  explanations,  the  interested  reader  is 
referred  to  references  24  and  25. 


2.1.1  MPIE-MoM 


To  calculate  the  radiating  and  induced  currents  of  the  RF  structure,  a  MoM  procedure  is  solved 
through  the  use  of  a  standard  form  of  the  electric  field  integral  equation.  For  an  arbitrary  field 
impinging  upon  a  perfectly  conducting  structure,  the  mixed  potential  integral  equation  (MPIE), 
in  terms  of  vector  potential  A  and  scalar  potential  o ,  is  given  by  the  following: 


nx 


[Einc(r)  +  Escat(f)\= 0,  re  structure  surface  s; 


Escat{r)  =  -jojA(r)  -  VO (r) ; 


(1) 

(2) 


where 

A(7)  =  \G\r,r')-Js(r')dr ; 
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(3) 


O(r)  =  — —  f  G°(r,  r} )V-JS (r ' )dr' . 

j<°i 


(4) 


The  dyadic  and  scalar  Green’s  functions  have  been  derived  (26-28)  for  the  case  when  the  region 
of  interest  is,  for  example,  the  upper  medium  of  a  half-space: 
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In  the  above  expressions,  the  quantities  denoted  by  subscripts  “fs”  and  “i”  correspond  to  free 
space  and  image  components,  respectively.  The  Sommerfeld  integrals  (I’s)  needed  to  complete 
equations  5  and  6  are  defined  in  references  24  and  25.  A  set  of  expressions  similar  to 
equations  5-10  can  be  derived  for  the  case  when  the  region  of  interest  is  the  lower  medium  of  the 
half-space,  or  when  the  propagation  background  is  a  multi-layered  (horizontally  stratified) 
medium. 


After  meshing  the  structure  of  interest  with  triangular  facets,  by  expanding  the  unknown  surface 
current  density  Js  as  a  summation  of  Rao-Wilton-Glisson  basis  functions  (29,  30)  and  then 

testing  equation  1  with  the  same  type  of  functions,  a  matrix  equation  is  formed,  enabling  the 
calculation  of  the  surface  current  distribution  throughout  the  structure.  A  procedure  similar  to  the 
one  shown  in  reference  29  can  be  used  for  setting  up  the  matrix  equation;  and,  with  the 
introduction  of  normalized  area  coordinates,  the  face-pair  combination  evaluation  can  be 
implemented  to  speed  up  the  impedance  matrix  filling  process  by  eliminating  redundant 
computational  routines  inherent  in  the  edge-pair  combination  evaluation.  For  a  voltage- source 
excited  structure,  a  delta-gap  model  is  assumed  for  the  excitation  vector  in  the  matrix  equation; 
in  this  model,  an  incident  field  only  exists  over  the  triangular  patches  spanning  the  gap  location. 
For  a  more  general  incident  field,  the  value  of  the  incident  field  is  assumed  to  be  constant  over 
each  patch  and  equal  to  its  value  at  the  patch’s  centroid. 

2.1.2  Approximation  for  Green’s  Functions 

Note  that  two  primary  difficulties  are  encountered  when  calculating  the  impedance  matrix 
elements  using  equations  5  and  6.  One  is  due  to  the  singularities  that  would  emerge  in  the 
expressions  g/s,  gi,  gxz,  and  gyz  for  certain  impedance  elements;  however,  singularity  extraction 
techniques  have  been  used  successfully  in  obtaining  closed-form  solutions  to  the  so-called  self¬ 
terms  (31-33).  Also,  a  quadrature  routine  can  be  employed  in  computing  the  elements  containing 
gxz  and  gyz,  even  in  the  presence  of  singularities.  The  second  difficulty  concerns  the  calculation  of 
the  impedance  matrix  elements  in  the  presence  of  the  Sommerfeld  integrals.  Implementing  an 
exact  routine  to  calculate  the  impedance  matrix  elements  would  necessitate  the  evaluation  of  a 
triple  integral.  To  expedite  the  matrix  filling  process,  different  common  schemes  are  available  in 
efficiently  calculating  the  Sommerfeld  integrals — either  by  speeding  up  the  convergence  rate  of 
the  integrals  or  approximating  the  integrals  with  closed-form  expressions.  Since  the  Sommerfeld 
integrals  are  slowly  convergent  (especially  for  small  (z  +  z’))  along  the  original  path  on  the  real- 
axis,  a  deformed  path  (e.g.,  steepest  descent  path)  can  be  exploited  to  achieve  faster 
convergence.  To  avoid  repeated  calculation  of  the  same  integrals,  another  class  of  methods  uses 
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a  pre-compute/tabulate/interpolate  routine  in  which  the  values  of  the  integrals  at  a  finite  number 
of  sampled  spatial  locations  are  calculated  in  advance  (using  convergence  acceleration 
techniques)  and  then  are  stored  in  a  database  to  serve  as  a  base  set  for  interpolation  by  the  MoM 
procedure.  In  this  study,  the  method  of  discrete  complex  images  is  employed  to  approximate  all 
the  Sommerfeld  integrals  with  closed-form  analytical  expressions;  thus,  by  eliminating  the 
innermost  of  the  triple  integral,  the  remaining  double  integral  can  be  carried  out  readily  using  a 
standard  quadrature  algorithm.  After  mapping  the  integrands  as  a  sum  of  N  exponentials  of  the 
form  6-aa  on  the  complex  koz-plane  using  the  generalized  pencil  of  function  method  (34- 

v=l 

37),  either  of  the  following  two  identities  is  applied  in  approximating  the  resulting  simplified 
Sommerfeld  integrals  in  canonical  forms: 


J 


e-jkoz(z+z') 


Jo(kpP)kpdkp 


-jknR: 

e  0  1 


(11) 


j  = 
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p(z  +  z')e 


-iKA 


R. 


3  +  3  jkaRi 


(12) 


Following  the  recipe  given  above,  the  five  Sommerfeld  integrals  of  interest  have  been 
approximated  with  an  error  of  less  than  1%  in  this  study.  In  a  particular  propagation  medium, 
each  set  of  complex  images  is  only  valid  in  a  particular  spatial  and  spectral  domain — in  other 
words,  each  set  affords  accurate  estimation  only  over  a  limited  range  of  p,  (z  +  z’),  and 
frequency.  (Formally,  unlike  that  of  equation  11,  the  right  side  of  equation  12  is  not  in  the  proper 
form  to  be  labeled  as  a  “complex  image”  contribution;  nonetheless,  the  name  “method  of  discrete 
images”  is  adopted  here  when  referring  to  this  technique.) 

In  examining  equations  5  and  6,  it  can  be  seen  that  for  a  structure  supporting  purely  horizontal 
currents  (e.g.,  distributed  filter),  only  two  Sommerfeld  integrals  are  needed  in  setting  up  the 
matrix  equation.  Furthermore,  for  a  general  three-dimensional  (3-D)  structure,  it  has  been 
observed  that  the  non-diagonal  components  (i.e.,  the  directionally  mixed  response  terms  g)  ,  G)  , 

G) ,  and  g)  )  of  the  dyadic  Green’s  function  have  only  minor  contributions  to  the  total  impedance 
matrix. 


2.2  Wideband  Lumped-element  Equivalent  Circuit  Synthesis 

RF  components  such  as  transmission  lines,  antennas,  and  distributed  filters  are  translated  into 
equivalent  circuits  for  subsequent  simulation  in  a  SPICE  solver.  The  multi-port  S-parameters  of 
the  structures  are  first  obtained  with  the  MoM  technique  and  then  transformed  into  F-parameters. 
The  frequency-domain  transfer  function — that  is,  in  this  case,  the  function  representing  the 
admittance  parameter,  YEm(s) — for  each  component  can  be  expanded  as  a  rational  approximation 
of  the  form 
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(13> 

/  s-Pi 

where  s  =  jco  is  the  complex  frequency;  r,  and  pt  are  the  residues  and  poles,  respectively;  and  d,  e 
are  constants.  Note  that  r;  and  p,  appear  as  either  pure  real  numbers  or  complex  conjugate  pairs. 
After  the  unknown  parameters  in  equation  13  are  estimated  with  the  two-stage  vector  fitting 
approach  outlined  by  Gustavsen  (38-41),  a  wideband  equivalent  circuit  network  for  the  RF 
component  is  synthesized  by  converting  each  term  within  equation  13  into  a  lumped-element 
circuit  using  the  procedure  prescribed  in  reference  42. 

2.3  Circuit  Simulation 

The  response  of  the  overall  equivalent  network  for  the  RF  structure — consisting  of  the  wideband 
lumped-element  sub-circuits  derived  in  section  2.2  and  the  extracted  nonlinear  components — is 
simulated  in  a  SPICE  solver.  Specifically,  a  time-domain  transient  analysis  is  carried  out  to 
obtain  the  output  harmonic  amplitudes  for  an  excitation  at  a  given  fundamental  frequency.  (Note 
that  a  harmonic-balance-based  circuit  simulator  can  also  be  used  to  realize  a  frequency-domain 
characterization  [43].)  A  first-order  approximation  is  employed  to  deduce  the  effective  stimulus 
needed  for  the  network.  The  re-radiated  power  at  each  harmonic  can  be  found  by  evaluating  the 
power  dissipated  by  the  radiation  resistances  within  the  equivalent  circuit.  The  harmonic  current 
densities  on  the  structure,  which  are  needed  for  the  field  pattern  calculations,  are  approximated 
by  the  linear  MoM  solver  (after  properly  terminating  the  structure  ports  with  equivalent 
impedances). 

2.4  Electromagnetic  Field  Propagator 

Once  the  current  distributions  have  been  obtained  with  the  MPIE-MoM,  the  radiation  properties 
of  the  RF  structure  can  be  characterized  by  applying  a  propagation  model  in  accordance  with  the 
geometry  and  physical  properties  of  the  propagation  medium.  For  a  current  source  located  over  a 
flat  ground  surface,  a  second-order  asymptotic  approximation  can  be  derived  from  an  exact 
solution  in  computing  the  radiated  fields;  the  formulations  necessary  for  this  step  can  be  deduced 
from  the  expressions  given  in  appendix  A  specialized  to  a  half-space  medium.  It  should  be 
mentioned  that  although  the  hybrid  MoM  and  asymptotic  propagator  approach  as  outlined  above 
has  been  described  only  for  currents  located  above  the  air/ground  interface,  an  equally  valid  and 
conceptually  similar  set  of  formulations  and  procedures  exists  for  the  subsurface  case.  The  only 
major  deviation  results  from  the  fact  that  when  the  radiators  are  submerged  in  a  lossy  medium, 
the  mode  of  propagation  taken  by  the  fields  may  be  in  the  form  of  lateral  waves  instead  of 
geometric-optics  contributions  and  Norton  waves  (44,  45);  hence,  the  asymptotic  propagator 
should  be  modified  accordingly. 

For  through-wall  environments,  the  radiation  from  sources  located  inside  a  room  is  calculated 
with  a  physical-optics  and  volumetric  polarization  current-based  approach  by  following  these 
steps: 
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1.  The  incident  waves  at  the  inner  wall  surfaces  are  determined  by  the  half-space  asymptotic 
method  discussed  above. 

2.  The  incident  waves — approximated  as  plane  waves — are  decomposed  into  transverse 
electric  (TE)  and  transverse  magnetic  (TM)  modes.  By  treating  the  walls  as  standard  two- 
layer  dielectric  structures,  the  modal  fields  inside  the  walls  are  then  estimated  with  the 
plane-wave  transmission  coefficients  of  the  media. 

3.  The  polarization  currents  for  replacing  the  walls  are  derived  from  the  volumetric 
equivalence  principle  as 

Jp(T)=-i<oe,(er-i)Ew(Ty.  (14) 

where  Ew(r)  is  the  electric  field  inside  the  wall  and  er  is  the  complex  relative  dielectric 
constant  of  the  wall  material. 

4.  Subsequently,  by  using  the  half-space  asymptotic  propagator,  the  signal  outside  the  room  is 
found  by  summing  the  radiation  from  the  wall  polarization  currents  and  the  direct  radiation 
due  to  the  RF  structure  in  the  absence  of  the  room: 

Kuj?)*  \Ghs{f,r)-Jp(r)dr  +  J (r , r ') - (r 'yr ' ,  (15) 

wall  RF  strucure 

where  Ghs(r,r')  is  the  half-space  dyadic  Green’s  function  from  appendix  A. 


The  advantages  of  the  outlined  semi-analytic  approach  are  that  it  is  computationally  faster  and 
more  efficient  than  an  exact  full-wave  method,  and  it  has  the  ability  to  incorporate  surface-wave 
contributions  into  the  solution.  The  shortcoming  is  that  the  technique  only  provides  an 
approximate  solution:  only  the  first-order  wall  reflection/scattering  effects  are  included  in  the 
solution.  (Note  that  multiple  reflections  within  the  walls  are  taken  into  account  for  the  calculation 
of  the  fields  in  step  2.)  To  improve  the  accuracy  of  the  solver,  an  iterative  physical-optics  routine 
can  be  used  to  include  shadowing  and  higher-order  effects. 


3.  Numerical  Examples:  Applications  of  the  Hybrid  Solver 


3.1  Circuit-Component-Loaded  Scatterer  in  Half-space  Environment 

As  most  RF  transceivers  contain  an  antenna  and  a  filter  element  at  their  front-ends,  the  structure 
in  figure  4  is  chosen  for  illustrative  purposes  in  this  section.  The  filter  is  a  microstrip  dual-band 
bandpass  filter  with  operational  bands  at  -900  MHz  and  -2000  MHz.  The  dual-band  topology 
results  from  the  use  of  the  combination  of  a  broadband  shunt  stub  bandpass  filter  and  an  open- 
stub  bandstop  filter.  The  broadband  bandpass  and  bandstop  filters  basically  composed  of 
approximately  quarter-wavelength  long  stubs  spaced  by  quarter-wavelength  transmission  line 
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sections.  Separately,  the  two  filters  can  be  translated  into  transmission  line-based  structures  by 
applying  Richard’s  transformation  to  their  respective  lumped-element-based,  low- 
pass/bandpass/bandstop  prototypes.  The  antenna — connected  to  the  microstrip  filter  and  printed 
on  the  other  side  of  the  two-layer  substrate — is  a  dual-band,  slot-loaded  microstrip  circular  patch, 
with  the  feeding  point  at  an  off-center  position.  A  shorting  connection  is  added  at  the  center  to 
reduce  the  size  of  the  patch.  Essentially,  the  patch  behaves  as  a  quarter-wave  radiator,  and  the 
fundamental  mode  of  the  circular  cavity  provides  the  resonance  at  the  lower  frequency  band, 
whereas  the  slot  supplies  the  resonance  at  the  higher  frequency  band.  (Note  that  dual-band 
performance  for  the  center-shorted  circular  cavity  can  also  be  established  without  the  slot. 
However,  it  is  seen  in  the  initial  design  that  the  second  resonance  frequency  of  the  cavity  is 
higher  than  predicted.  The  presence  of  the  slot  enables  the  generation  of  the  needed  response  at 
-1900-2000  MHz  while  maintaining  a  compact  geometry.  Consequently,  the  second  resonance 
of  the  cavity  is  affected  by  the  slot,  but  the  first  resonance  is  not.)  To  introduce  a  nonlinear 
element  into  the  structure,  a  varactor  is  added  to  the  antenna,  providing  capacitive  tuning 
capability.  The  effective  capacitance  of  the  varactor  as  a  function  of  the  DC  bias  voltage  is 
shown  in  figure  5.  Simulations  of  the  filter  and  the  antenna  are  carried  out  using  the  linear  MoM 
approach  from  section  2 — assuming  infinite,  layered  substrates.  The  S-parameters  of  the 
components  are  shown  in  figures  6  and  7. 
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Figure  4.  Antenna  and  filter  for  RF  front-end:  (a)  top  view;  (b)  bottom  view; 
and  (c)  side  view.  sri  =  sr2  =  4.5  (tan8  ~  0.04);  the  thickness  of 
each  substrate  is  1.6  mm.  The  distributed  filter  has  dimensions 
-10.7  cm  x  7  cm;  the  radius  of  the  circular  patch  antenna  is  1.79  cm, 
and  the  length  of  the  slot  is  1.14  cm. 
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DC  Bias  Voltage  [V] 


Figure  5.  Effective  capacitance  of  varactor  used  in  the  setup  shown  in 
figure  4  as  function  of  DC  bias  voltage. 


Frequency  [MHz] 


Figure  6.  \Sn\  for  the  dual-band,  slot-loaded  circular  patch  antenna  as 
function  of  frequency  and  varactor  capacitance. 
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Figure  7.  ^-parameters  for  microstrip  dual-band  bandpass  filter  as 
function  of  frequency. 

The  scattering  from  the  RF  structure  in  figure  4  is  calculated  with  the  hybrid  approach  put  forth 
in  section  2.  Note  that  for  the  analysis  that  follows,  it  is  assumed  that  the  filter  input  has  been 
terminated  with  a  50  Q  load  and  the  varactor  is  biased  at  0  V.  The  structure  is  placed  at  10  cm 
below  the  ground  surface  (the  ground  has  relative  dielectric  constant  er’  =  4  and  conductivity 
era  =  10  mS/m).  The  excitation  is  provided  by  a  vertical  electric  dipole — located  10  m  from  the 
scatterer  at  a  2-m  height — transmitting  a  monochromatic  wave  with/0  =  840  MHz.  Figure  8 
shows  the  backscattered  responses  of  the  scatterer  at  the  fundamental  and  the  second  harmonic 
for  two  levels  of  transmitted  power — that  is,  for  Ptx  =  100  W  and  1000  W.  The  received  power 
at  the  receiver  is  estimated  by  first  propagating  the  fields  from  the  scatterer  using  the  asymptotic 
solver  and  then  calculating  the  power  intercepted  by  a  polarization-matched  receiving  antenna 
with  0  dB  gain.  The  near-quadratic  dependence  of  the  scattered  power  at  2/0  on  the  transmitted 
power  is  evident.  The  instability  of  the  time-domain  SPICE  solver  prevents  the  responses  at  the 
higher-order  harmonics  to  be  derived  here;  however,  note  that  an  accurate  spectrum 
characterization  may  be  found  using  a  harmonic-balance-based  circuit  simulator.  Figure  9  shows 
the  variation  of  the  scattered  fields  as  a  function  of  the  azimuth  angle. 
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Figure  8.  Backscattered  response  of  scatterer  in  figure  4  at  fundamental  and  second  harmonic 
for  two  levels  of  transmitted  power.  The  noise  floor  ( kT)  shown  is  at  -173.8  dBm/Hz. 


Figure  9.  Variation  of  scattered  electric  field  components  as 
a  function  of  azimuth  angle. 

To  investigate  the  decay  of  the  signal  as  a  function  of  sensing  distance,  the  backscattered 
responses  are  calculated  for  three  different  standoff  ranges — r  =  10  m,  30  m,  and  100  m.  The 
results  in  figure  10  demonstrate  that  the  fundamental  scattering  component  falls  off  as  oc  r“8 , 
whereas  the  second  harmonic  falls  off  as  oc  r  12  .  These  behaviors  are  consistent  with  theoretical 
results  derived  using  analytical  models  for  the  nonlinear  circuit  response  and  field  propagation. 
As  a  generalization,  in  the  asymptotic  regime,  it  can  be  shown  that — for  the  backscattering-based 
sensing  scenario  in  which  the  interaction  points  are  positioned  close  to  the  ground  surface — the 
«-th  harmonic  return  approximately  depends  on  the  distance  as  oc  r  4("+l1 . 
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Figure  10.  Backscattered  response  as  a  function  of  range:  (a )f0  and  (b)  2 f„. 

To  see  the  effects  of  a  more  complicated  nonlinear  load  on  the  scattering  response,  a  low-noise, 
traveling-wave  amplifier  is  connected  at  the  filter  input  port.  The  design  and  circuit  layout  for  the 
amplifier  are  considered  in  detail  within  appendix  C;  a  two-stage  topology  (figure  C-6a)  is 
employed  for  the  scattering  analysis  here.  The  received  power  level  as  a  function  of  frequency — 
up  to  the  third  harmonic — is  displayed  in  figure  1 1  for  the  monochromatic  incidence  case;  the 
simulation  parameters  are  the  same  as  those  of  figure  8.  It  should  be  mentioned  that,  for  this 
example  problem,  in  addition  to  the  amplifier,  more  components — passive  and  active,  linear  and 
nonlinear — can  be  connected  to  the  RF  structure.  These  components  are  “transparent”  to  the 
electromagnetic  solver — that  is,  they  are  replaced  only  by  equivalent  loads  in  the  field  analysis, 
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for  calculating  the  effective  system  stimulus  and  the  structure  current  densities;  full  circuit 
models  of  these  components,  however,  can  be  included  in  the  SPICE  solver  to  characterize  the 
circuit  spectral  response. 


Figure  1 1 .  Backscattered  responses  at  fundamental,  second,  and  third  harmonics 
of  scatterer  in  figure  4  connected  to  a  traveling-wave  amplifier. 

Once  the  scattered  fields  from  the  RF  structure  are  obtained,  an  imaging  algorithm  can  be  used  to 
process  the  return  signals  for  localization.  In  this  work,  a  time-reversal-based  MUSIC  (TR- 
MUSIC)  algorithm  is  applied  to  generate  the  imaging  functional.  Consider  the  multi-static 
sensing  of  a  scene  with  M  nonlinear  scatterers  using  an  iV-element  array,  the  scattered  field 
response  (at  frequency  oos)  received  at  the  array  due  to  excitation  (at  frequency  coi)  by  the  n-th 
transmitter  of  the  array  can  be  shown  to  be 
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where  rn  (n=  1,2,  N)  is  the  location  of  the  n-th  array  element,  rsm  (m  =  1,2,  . . .,  M)  is  the 
location  of  the  m-th  scatterer,  <ysm  is  the  scattering  strength  of  the  m-th  scatterer,  and 

G(r,r\  (o)  is  the  Green’s  function  of  the  environment.  From  equation  16,  it  is  seen  that  the 
scattered  field  signal  subspace  is  spanned  by  the  Green’s  function  vectors  G(rsl,cos),  G ( r 2 , &>s ) , 
. . .,  G(rsM, cos ) — in  other  words,  the  scatterer  locations  are  encoded  within  the  subspace 
representation  of  the  received  signal  (46).  In  more  concise  form,  after  invoking  the  reciprocity 
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principle,  the  frequency-domain  NxN  multi-static  scattered  signal  matrix  for  the  array  can  be 
written  as 


M 

K {o)s,o)i)G(fsm,o)s)GT  (rim,rq.) ,  (17) 

m=l 

in  which  T  represents  the  transpose  operation,  and  the  matrix  element  Kpq(o)s,  cy)  is  the  response 
at  the  p-th  array  element  due  to  excitation  at  the  g-th  array  element.  In  practical  applications, 
K(a>s,  CO;)  is  simply  taken  as  the  measurement  matrix.  For  image  generation,  the  signal  subspace 
of  K{ (os,  o)j)  is  computed  using  a  singular  value  decomposition  procedure: 

K  (cos ,  cot )  =  U  (o)s ,  <»,. )  A (a>s ,  ) V  (cos ,  cot  )H ,  (18) 


where  U( G)s,  Wi)  and  V(cos,  ay)  are  unitary  matrices,  A  (ay,  oy)  contains  the  singular  values  of 
K{ (os,  Wi),  and  H  denotes  the  conjugate  transpose  operation.  Noting  that  the  column  vectors  of 
£/(oy,  coi)  supply  the  eigenvectors  up  ( cos ,  01 )  (p  =  1,2,  ...,  AO,  and  assuming  that  the  scattered 

signal  subspace  is  spanned  by  the  eigenvectors  corresponding  to  the  first  L  significant  (or  non¬ 
zero)  eigenvalues  and  the  null  subspace  is  spanned  by  the  remaining  eigenvectors  corresponding 
to  the  non-significant  (or  zero)  eigenvalues,  an  imaging  functional  can  be  constructed  as 


(  N  _ 

24 
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J 

where  the  Green’s  function  vector  G(f,a>s )  can  be  computed  using  numerical  or  analytical 
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methods.  The  imaging  functional  in  equation  19  peaks  at  the  scatterer  locations.  In  effect,  then, 
an  image  of  the  scene  can  be  derived  using  the  scattering  response  at  a  single  receiving 
frequency.  To  increase  the  stability  of  the  algorithm,  equation  19  can  be  performed  over  multiple 
receiving  frequencies  (e.g.,  harmonics)  or  a  band  of  receiving  frequencies: 


(  N 

2  4 

Q(r,  Q,ru,)  = 

M 

M 

(up(o)s,o)i),G(r,cos)) 

^  Q  p=L[cos  ,0)1)+! 

J 

(20) 


in  which  Q  denotes  the  range  of  available  receiving  frequencies.  Equations  19  and  20  can  also  be 
generalized  to  accommodate  scenarios  in  which  multiple  excitation  frequencies  (multi-coi)  are 
utilized.  (Note  that  TR-MUSIC  based  on  the  conventional,  or  linear,  scattering  response  is 
derived  by  setting  ws  =  cy  in  equations  16-20.) 

The  imaging  functional  in  equation  19  is  employed  for  imaging  the  scenario  displayed  in 
figure  12.  The  scene  consists  of  the  circuit-component-loaded  scatterer  from  figure  4  buried  in 
the  ground,  along  with  two  (linear)  clutter  objects.  The  sensing  array  is  composed  of  N  =  16 
transceivers  distributed  over  a  2-m  wide  aperture  and  the  standoff  range  is  ~6.6  m.  The  harmonic 
multi-static  scattered  signal  matrices  of  the  scene  at  f0  and  2/0  are  calculated  with  the  hybrid 


17 


solver.  G(r,  <ds  )  in  equation  19  is  calculated  with  the  formulations  in  appendix  A  (or  reference 

24).  The  images  at  the  two  frequencies  are  included  in  figure  12.  Note  that  the  target  is 
accurately  localized  in  both  images,  but  the  clutter  signal  components  do  not  show  up  in  the 
image  at  2/0. 


Figure  12.  Imaging  results  for  a  buried  target  scenario:  (a)  Scene  containing  target  and  clutter,  (b)  image 

generated  using  scattering  response  at  f0  =  840  MHz  (both  the  target  and  clutter  are  imaged),  and 
(c)  image  generated  using  scattering  response  at  2/0  =  1680  MHz  (only  the  target  is  captured). 

3.2  Circuit-Component-Loaded  Scatterer  in  Room  Environment 

To  investigate  the  response  of  a  circuit-component-loaded  scatterer  for  through-wall  sensing 
applications,  the  RF  structure  from  section  3.1  is  placed  in  the  room  environment  shown  in 
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figure  2.  In  this  example,  the  walls,  ground,  and  roof  (not  shown  in  figure  2)  are  homogenous 
media  with  er’  =  6  and  era  =  30  mS/m.  Note  that  as  before,  the  filter  input  has  been  terminated 
with  a  50  Q  load  and  the  varactor  is  biased  at  0  V.  The  RF  structure  is  placed  at  1  m  above  the 
ground;  the  excitation  is  provided  by  a  vertical  electric  dipole  located  outside  the  room  at  a 
height  of  1  m,  transmitting  at  f0  =  840  MHz,  with  Pjx  =  100  W;  and  a  receiving  antenna  (co¬ 
located  with  the  transmitter)  with  0  dB  gain  is  assumed.  The  hybrid  solver  outlined  in  section  2  is 
applied  for  deriving  the  backscattered  responses  at/0  and  2 f0.  To  calculate  the  initial  incident 
field  on  the  RF  structure,  the  reciprocity  principle  is  exploited:  the  locations  of  the  transmitter 
and  observation  point  are  interchanged  and  then  the  fields  are  propagated  with  the  physical- 
optics  and  volumetric  polarization  current-based  approach  from  section  2.4.  Figure  13  shows  the 
backscattered  responses  of  the  scatterer  at  the  fundamental  and  the  second  harmonic.  The  bistatic 
scattering  pattern  as  a  function  of  azimuth  angle  is  shown  in  figure  14. 
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Figure  13.  Backscattered  response  of  scenario  shown  in  figure  2:  (a)  without 
and  (b)  with  circuit-component-loaded  scatterer  in  room. 
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Figure  14.  Backscattered  response  at  fundamental  and  second 
harmonic  as  function  of  azimuth  angle  for  scenario 
in  figure  2. 

In  order  to  verify  the  accuracy  of  the  field  propagator  in  the  presence  of  a  wall,  simulation  results 
are  compared  with  measurement  data  obtained  by  Dagefu  and  Sarabandi  (47).  The  wall  is  solid 
concrete  with  eT’  =  4  and  rrd  =  20  mS/m  and  has  dimensions  2.56  m  x  1.56  m  x  0.09  m.  Two 
setups  are  considered.  In  the  first  case  (figure  15a),  the  transmitter  is  located  1.74  m  from  the 
wall  at  a  height  of  0.87  m;  the  signal  strength  along  an  observation  line  0.91  m  from  the  wall  at 
0.79  m  height  is  recorded.  In  the  second  case  (figure  16a),  the  transmitter  is  located  1.74  m  from 
the  wall  at  a  height  of  0.24  m;  the  signal  strength  along  an  observation  line  perpendicular  to  the 
wall  at  0.20  m  height  is  recorded.  In  this  example,  for  the  electromagnetic  simulations,  f0  = 

2  GHz.  Note  that  the  measurement  scenario  only  considers  the  linear  propagation  effects  at  the 
fundamental  frequency  of  interest.  Figures  15b  and  16b  show  that  there  is  good  agreement 
between  the  analytically  derived  results  and  the  measured  data. 
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Figure  15.  Simulated  vs.  measured  responses  for  propagation  in  presence  of  a  wall. 
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Figure  16.  Simulated  vs.  measured  responses  for  propagation  in  presence  of  a  wall. 


4.  Conclusions 


To  study  the  radio  harmonic  scattering  responses  of  circuit-component-loaded  objects,  a  hybrid 
electromagnetic/circuit  solver  is  developed  in  this  work.  The  proposed  numerical  analysis 
approach  exploits  a  combination  of  3-D  electromagnetic  simulations  (based  on  full- wave  and 
asymptotic  techniques)  in  the  propagation  and  RF  component  domains  and  SPICE  simulations  in 
the  circuit  domain.  In  sum,  a  systematic  strategy — one  which  is  capable  of  tracking  the  signal 
interactions  from  the  transmitter  terminal  to  the  receiver  terminal — is  established  for  attaining 
generalized  solutions  of  the  multiscale  harmonic  scattering  problem.  This  terminal-to-terminal 
(or  end-to-end)  solver  has  been  applied  to  two  common  radar  sensing  problems  of  interest.  In 
particular,  considered  herein  in  detail  are  the  linear  and  nonlinear  scattering  effects  from  an  RF 
circuit-loaded  object  placed  in  either  near-ground  or  through-wall  settings.  The  scattering 
responses  as  functions  of  excitation  power,  sensing  range,  and  azimuth  angle  are  examined,  and 
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the  localization,  or  imaging,  of  the  scatterer  using  both  the  linear  and  the  nonlinear  returns  has 
also  been  investigated.  Through  the  theoretical  developments  and  illustrative  examples  put  forth, 
it  is  demonstrated  that  the  proposed  hybrid  solver — though  not  completely  full-wave-based — 
establishes  a  practical  and  computationally  efficient  solution  process  for  characterizing  harmonic 
scattering  in  complex  environments. 


24 


5.  References 


1.  Piket-May,  M.;  Taflove,  A.;  Baron,  J.  FD-TD  Modeling  of  Digital  Signal  Propagation  in  3-D 
Circuits  With  Passive  and  Active  Loads.  IEEE  Trans.  Microw.  Theory  and  Techniques  Aug. 
1994,  43  (8),  1514-1523. 

2.  Sui,  W.;  Christensen,  D.  A.;  Durney,  C.  H.  Extending  the  Two-Dimensional  FD-TD  Method 
to  Hybrid  Electromagnetic  Systems  With  Active  and  Passive  Lumped  Elements.  IEEE 
Trans.  Microw.  Theory  and  Techniques  Apr.  1992,  40,  726-730. 

3.  Ciampolini,  P.;  Mezzanotte,  P.;  Roselli,  L.;  Sorrentino,  R.  Accurate  and  Efficient  Circuit 
Simulation  With  Lumped-Element  FDTD  Technique.  IEEE  Trans.  Microw.  Theory  and 
Techniques  Dec.  1996,  44  (12),  2207-2215. 

4.  Pereda,  J.  A.;  Alimenti,  F.;  Mezzanotte,  P.;  Roselli,  L.;  Sorrentino,  R.  A  New  Algorithm  for 
the  Incorporation  of  Arbitrary  Linear  Lumped  Networks  into  FDTD  Simulators.  IEEE 
Trans.  Microw.  Theory  and  Techniques  Jun.  1999,  47  (6),  943-949. 

5.  Huang,  C.-C.;  Chu,  T.-H.  Analysis  of  Wire  Scatterers  with  Nonlinear  or  Time-Harmonic 
Loads  in  the  Frequency  Domain.  IEEE  Trans.  Antennas  and  Propagation  Jan.  1993,  41  (1), 
25-30. 

6.  Parron,  J.;  Collado,  C.;  Mateu,  J.;  Rius,  J.  M.;  Duffo,  N.;  O'Callaghan,  J.  M.  General 
Electromagnetic  Simulation  Tool  to  Predict  the  Microwave  Nonlinear  Response  of  Planar, 
Arbitrarily-Shaped  HTS  Structures.  IEEE  Trans.  Applied  Superconductivity  Mar.  2001, 11 
(1),  399-102. 

7.  Liu,  T.  K.;  Tesche,  F.  M.  Analysis  of  Antennas  and  Scatterers  with  Nonlinear  Loads.  IEEE 
Trans.  Antennas  and  Propagation  Mar.  1976,  24,  131-139. 

8.  Liu,  T.  K.;  Tesche,  F.  M.;  Deadrick,  F.  J.  Transient  Excitation  of  an  Antenna  with  a 
Nonlinear  Load:  Numerical  and  Experimental  Results.  IEEE  Trans.  Antennas  and 
Propagation  Jul.  1977,  25,  539-542. 

9.  Sarkar,  T.  K.;  Weiner,  D.  D.  Scattering  Analysis  of  Nonlinearly  Loaded  Antennas.  IEEE 
Trans.  Antennas  and  Propagation  Mar.  1976,  24,  125-131. 

10.  Zaboronkova,  T.  M.;  Gorbachev,  A.  A.;  Vasenkov,  A.  A.  Nonlinear  Electromagnetic  Wave 
Scattering  by  Loaded  Dipole  Antennas.  Applied  Electromagnetism,  Proc.  Second 
International  Symposium  of  Trans.  Black  Sea  Region,  pp.  5-6,  2000. 

11.  Kanda,  M.  Analytical  and  Numerical  Techniques  for  Analyzing  an  Electrically  Short  Dipole 
With  a  Nonlinear  Load.  IEEE  Trans.  Antennas  and  Propagation  Jan.  1980,  28,  71-78. 


25 


12.  Luebbers,  R.;  Beggs,  J.;  Chamberlin,  K.  Finite  Difference  Time-Domain  Calculation  of 
Transients  in  Antennas  with  Nonlinear  Loads.  IEEE  Trans.  Antennas  and  Propagation  May 
1993,  41  (5),  566-573. 

13.  Porti,  J.  A.;  Morente,  J.  A.  A  Numerical  Analysis  of  Wire  Antennas  Loaded  With  Varistor- 
Composite  Materials.  IEEE  Trans.  Electromagn.  Compat.  Feb.  1994,  36  (1),  23-31. 

14.  Lee,  K.  C.  An  Efficient  Algorithm  for  the  Steady-State  Analysis  of  a  Nonlinearly  Loaded 
Antenna  Array.  J.  Electromagn.  Waves  Applicat.  2000, 14,  1373-1382. 

15.  Analysis  of  Large  Nonlinearly  Loaded  Antenna  Arrays  under  Multi  tones  Excitation,” 
Microw.  Opt.  Technol.  Lett.  Jun.  2000,  319-323. 

16.  Two  Efficient  Algorithms  for  the  Analyses  of  a  Nonlinearly  Loaded  Antenna  and  Antenna 
Array  in  the  Frequency  Domain.  IEEE  Trans.  Electromagn.  Compat.  Nov.  2000,  45  (4), 
339-346. 

17.  Genetic  Algorithms  Based  Analyses  of  Nonlinearly  Loaded  Antenna  Arrays  Including 
Mutual  Coupling  Effects,”  IEEE  Trans.  Antennas  and  Propagation  Apr.  2003,  51  (4),  776- 
781,. 

18.  Kalialakis,  C.;  Hall,  P.  S.  Analysis  and  Experiment  on  Harmonic  Radiation  and  Frequency 
Tuning  of  Varactor-Loaded  Microstrip  Antennas.  IET  Microw.  Antennas  Propagation  Apr. 
2007, 1  (2),  527-535. 

19.  Zheng,  Y.;  Hristov,  A.;  Giere,  A.;  Scheele,  P.;  Jakoby,  R.  Analysis  on  Harmonic  Radiation 
of  Capacitively  Loaded  Tunable  PIFA  Based  on  Equivalent  Circuit.  Frequenz,  2007,  61, 
213-216. 

20.  Cui,  J.;  Zhang,  M.  Transient  Analysis  of  PCB  EMC  Problem.  8th  International  Symposium 
on  Antennas  Propagation  and  EM  Theory,  ISAPE  2008,  pp.  Ill  1-1 1 14,  Nov.  2008. 

21.  Kung,  F.;  Chuah,  H.-T.  A  Finite-Difference  Time-Domain  (FDTD)  Software  for  Simulation 
of  Printed  Circuit  Board  (PCB)  Assembly.  Progress  In  Electromagnetics  Research  2005,  50, 
299-335. 

22.  £akir,  G.;  £akir,  M.;  Sevgi,  L.  Electromagnetic  Radiation  from  Multilayer  Printed  Circuit 
Boards:  A  3D  FDTD-Based  Virtual  Emission  Predictor.  Turk  J  Elec  Eng  and  Comp  Sci 
2009,  77(3),  315-325. 

23.  Asai,  H.;  Watanabe,  T.;  Sasaki,  T.;  Araki,  K.  An  EMI  Simulator  Based  on  the  Parallel- 
Distributed  FDTD  Method  for  Large-Scale  Printed  Wiring  Boards.  Proc.  6th  IEEE 
Workshop  on  Signal  Propagation  on  Interconnects,  pp.  141-144,  May  2002. 

24.  Liao,  D.  H.  Physics-Based  Near-Earth  Radiowave  Propagation  Modeling  and  Simulation, 
Ph.D.  Thesis,  The  University  of  Michigan,  2009. 


26 


25.  Liao,  D.  H.;  Sarabandi,  K.  Terminal-to-Terminal  Hybrid  Full-Wave  Simulation  of  Low- 
Profile,  Electrically-Small,  Near-Ground  Antennas.  IEEE  Trans.  Antennas  and  Propagation 
Mar.  2008,  56  (3),  806-814. 

26.  Michalski,  K.  A.;  Zheng,  D.  Electromagnetic  Scattering  and  Radiation  by  Surfaces  of 
Arbitrary  Shape  in  Layered  Media,  Part  I:  Theory.  IEEE  Trans.  Antennas  and  Propagation 
Mar.  1990,  38,  335-344. 

27.  Electromagnetic  Scattering  and  Radiation  by  Surfaces  of  Arbitrary  Shape  in  Layered  Media, 
Part  II:  Implementation  and  Results  for  Contiguous  Half-Spaces.  IEEE  Trans.  Antennas  and 
Propagation  Mar.  1990,  38,  345-352. 

28.  Michalski,  K.  A.  On  the  Scalar  Potential  of  a  Point  Charge  Associated  with  a  Time- 
Harmonic  Dipole  in  a  Layered  Medium.  IEEE  Trans.  Antennas  and  Propagation  Nov.  1987, 
35,  1299-1301. 

29.  Rao,  S.  M.  Electromagnetic  Scattering  and  Radiation  of  Arbitrarily-Shaped  Surfaces  by 
Triangular  Patch  Modeling,”  Ph.  D.  dissertation,  The  University  of  Mississippi,  Aug.  1980. 

30.  Rao,  S.  M.;  Wilton,  D.  R.;  Glisson,  A.  W.  Electromagnetic  Scattering  by  Surfaces  of 
Arbitrary  Shape.  IEEE  Trans.  Antennas  and  Propagation  May  1982,  30,  409^-18. 

31.  Eibert,  T.  F.;  Hansen,  V.  On  the  Calculation  of  Potential  Integrals  for  Linear  Source 
Distributions  on  Triangular  Domains.  IEEE  Trans.  Antennas  and  Propagation  Dec.  1995, 

43,  1499-1502. 

32.  Graglia,  R.  D.  On  the  Numerical  Integration  of  the  Linear  Shape  Functions  Times  the  3-D 
Green’s  Function  or  Its  Gradient  on  a  Plane  Triangle.  IEEE  Trans.  Antennas  and 
Propagation  Oct.  1993,  41,  1448-1455. 

33.  Rao,  S.  M.;  Glisson,  A.  W.;  Wilton,  D.  R.;  Vidula,  B.  S.  A  Simple  Numerical  Solution 
Procedure  for  Statics  Problems  Involving  Arbitrary-Shaped  Surfaces.  IEEE  Trans.  Antennas 
and  Propagation  Dec.  1979,  27,  604-608. 

34.  Hua,  Y.;  Sarkar,  T.  K.  Generalized  Pencil-of-Function  Method  for  Extracting  Poles  of  an 
EM  system  from  Its  Transient  Response.  IEEE  Trans.  Antennas  and  Propagation  Feb.  1989, 
37,  229-234. 

35.  Sarkar,  T.  K.;  Pereira,  O.  Using  the  Matrix  Pencil  Method  to  Estimate  the  Parameters  of  a 
Sum  of  Complex  Exponentials,”  IEEE  Antennas  and  Propagation  Magazine  Feb.  1995,  37, 
48-55. 

36.  Shuley,  N.  V.;  Boix,  R.  R.;  Medina,  F.;  Homo,  M.  On  the  Fast  Approximation  of  Green’s 
Functions  in  MPIE  Formulations  for  Planar  Layered  Media.  IEEE  Trans.  Microw.  Theory 
and  Techniques  Sep.  2002,  50,  2185-2192. 


27 


37.  Hellen,  M.  K.;  Craddock,  I.  J.  Improving  the  Accuracy  of  the  Complex  Image  Technique. 
Microw.  Optical  Technology  Letters  Mar.  2001,  28,  402-406. 

38.  Gustavsen,  B.;  Semiyen,  A.  Rational  Approximation  of  Frequency  Domain  Responses  by 
Vector  Fitting.  IEEE  Trans.  Power  Delivery  Jul.  1999, 14  (3),  1052-1061. 

39.  Gustavsen,  B.  Improving  the  Pole  Relocating  Properties  of  Vector  Fitting.  IEEE  Trans. 
Power  Delivery  Jul.  2006,  21  (3),  1587-1592. 

40.  Morched,  A.;  Gustavsen,  B.;  Tartibi,  M.  A  Universal  Model  for  Accurate  Calculation  of 
Electromagnetic  Transients  on  Overhead  Lines  and  Underground  Cables.  IEEE  Trans. 
Power  Delivery  Jul.  1999,  14  (3),  1032-1038. 

41.  Gustavsen,  B.  Wide  Band  Modeling  of  Power  Transformers.  IEEE  Trans.  Power  Delivery 
Jan.  2004, 19  (1),  414^122. 

42.  Antonini,  G.  SPICE  Equivalent  Circuits  of  Frequency-Domain  Responses.  IEEE  Trans. 
Electromagnetic  Compatibility  Aug.  2003,  45  (3),  502-512. 

43.  Kundert,  K.  S.;  Sangiovanni-Vincentelli,  A.  Simulation  of  Nonlinear  Circuits  in  the 
Frequency  Domain.  IEEE  Trans.  Computer-Aided  Design  of  Integrated  Circuits  and 
Systems  Oct.  1986,  5  (4),  521-535. 

44.  Liao,  D.  H.;  Sarabandi,  K.  Near-Earth  Wave  Propagation  Characteristics  of  Electric  Dipole 
in  Presence  of  Vegetation  or  Snow  Layer.  IEEE  Trans.  Antennas  and  Propagation  Nov. 
2005,  53,  3747-3756. 

45.  Sarabandi,  K.;  Koh,  I.  Effect  of  Canopy-Air  Interface  Roughness  on  HF-VHF  Wave 
Propagation  in  Forest.  IEEE  Trans.  Antennas  and  Propagation.  Feb.  2002,  50,  111-121. 

46.  Lev-Ari,  H.;  Devaney,  A.  J.  The  Time-Reversal  Technique  Re-Interpreted:  Subspace-Based 
Signal  Processing  For  Multi-Static  Target  Location.  Proc.  IEEE  Sensor  Array  and 
Multichannel  Signal  Processing  Workshop,  pp.  509-513,  2000. 

47.  Dagefu,  F.  T.;  Sarabandi,  K.  Analysis  and  Modeling  of  Near-Ground  Wave  Propagation  in 
the  Presence  of  Building  Walls.  IEEE  Trans.  Antennas  and  Propagation  Jun.  2011,  59  (6), 
2368-2378. 

48.  Felsen,  L.  B.;  Marcuvitz,  N.  Radiation  and  Scattering  of  Waves',  NJ:  Prentice-Hall,  1973. 

49.  Jordan,  E.  C.;  Balmain,  K.  G.  Electromagnetic  Waves  and  Radiating  Systems',  Englewood 
Cliffs,  NJ:  Prentice-Hall,  1968. 

50.  Brekhovskikh,  L.  M.  Waves  in  Layered  Media',  NY:  Academic  Press,  1960. 

51.  Tsang,  L.;  Kong,  J.  A.;  Shin,  R.  T.  Theory  of  Microwave  Remote  Sensing;  NY:  Wiley,  1985. 


28 


52.  Kong,  J.  A.;  Tsang,  L.;  Simmons,  G.  Geophysical  Surface  Probing  With  Radio-Frequency 
Interferometry.  IEEE  Trans.  Antennas  and  Propagation.  Nov.  1974,  22. 

53.  Chew,  W.  C.  Waves  and  Fields  in  Inhomogeneous  Media',  NY:  Van  Nostrand  Reinhold, 
1990. 

54.  Beyer,  J.  B.;  Prasad,  S.  N.;  Becker,  R.  C.;  Nordman,  J.  E.;  Hohenwarter,  G.  K.  MESFET 
Distributed  Amplifier  Design  Guidelines.  IEEE  Trans.  Microw.  Theory  and  Techniques 

Mar.  1984,  32. 

55.  D.  M.  Pozar,  Microwave  Engineering',  2nd  Ed.,  John  Wiley  and  Sons,  NY,  1997. 


29 


Intentionally  Left  Blank. 


30 


Appendix  A.  Higher-Order  Asymptotic  Evaluation  of  Near-Earth  Wave 
Propagation  Characteristics  of  Electric  Current  Sources 


Analytical  formulations  are  derived  for  calculating  the  electric  fields  of  a  short  dipole  of  arbitrary 
orientation  situated  above  a  finitely  conducting  ground  plane  in  the  presence  of  a  dielectric  layer. 
The  discussion  presented  below  is  originally  from  the  work  by  Liao  (24).  Note  that  specialized 
results  for  the  half-space  medium  can  be  reached  by,  for  example,  setting  the  dielectric  properties 
of  the  dielectric  layer  to  those  of  air  (or  free  space). 

There  are  three  primary  configurations  of  interest,  as  a  result  of  different  combinations  of  the 
locations  of  the  electric  dipole  and  observation  point.  The  dipole  is  assumed  to  be  located  on  the 
z-axis  with  current  moment  vector  ij=io(i Xx  +  lyy  +  l,z)  >  and  the  observation  point — identified 

with  the  usual  cylindrical  coordinates — is  at  (p,  <j>,  z.).  The  vy-plane  is  aligned  with  the  top  of  the 
dielectric  layer,  itself  has  thickness  d.  In  the  type  of  propagation  problems  considered  here,  for 
analysis  restricted  to  the  far  field,  it  can  be  assumed  that  the  magnitudes  of  z\  z,  and  d  are  much 
smaller  than  that  of  p — the  radial  distance  between  the  dipole  and  the  observation  point.  Also 
note  that  primed  quantities  are  associated  with  the  source,  and  the  exp (-icot)  convention  is 
assumed  and  suppressed  in  all  the  formulations. 

A-l.  Case  1  (z’  >  0  and  z  >  0) 

Consider  the  simple  case  when  the  source  and  observation  points  are  both  located  in  the  upper 
layer  (air)  of  a  one-layer  medium,  the  total  electric  field  for  a  vertical  dipole  can  be  written  as 


tt1  total  _  direct  .  jp  scattered 

^  x,y,z  ^  x,y,z  ^  x,y,z 


(A-l) 


Specifically,  at  grazing  angle  (z  ’  «p  and  z  «  p),  each  component  of  the  field  can  be  expanded 
as  follows: 


E,olal  1  |  Rtm  ,  C,1  ,  c,2 

4  n  \p  p  p 2  p3 
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icoujle 


(z'-z)  (z'+z)Rtm  ,  cxj  Cx2 


2  2  3 
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,  i0MoIolzeik‘p  \{z'-z)  (z'+z)Rm  Cu  Cy-2 
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(A-2) 


(A-3) 


(A-4) 


The  first  term  in  the  brackets  is  the  direct  field;  the  second  term  is  the  geometrical-optics 
reflection  field;  and  the  rest  of  the  terms  account  for  higher-order  scattered  fields — where  the  Cs 
denote  unknown  factors.  Since  incidence  is  near  grazing  angle,  it  is  easy  to  show  that  Rtm  ~  -1, 
(1  -Rtm)  ~  2,  whereas  (1  -I -Rtm)  is  a  small  quantity  proportional  to  Up.  For  the  z-component  of  the 
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field,  it  can  be  seen  that  the  cancellation  of  the  direct  and  geometrical-optics  reflection  term 
generates  a  Up1  resultant — which  is  of  the  same  order  as  the  next  term  in  the  expansion. 
Likewise,  for  the  x  and  y  components,  the  sum  of  the  direct  and  geometrical-optics  reflection 
also  produces  a  term  that  is  of  the  same  order  as  the  next  higher-order  term.  Similar  analysis  can 
be  applied  to  a  horizontal  dipole.  Thus,  when  the  source  and  observation  points  are  close  to  the 
interface  (incidence  at  grazing  angle),  it  is  apparent  that  an  accurate  description  of  the  field 
components  necessitates  the  derivation  of  the  higher-order  term  that  comes  after  the  geometrical- 
optics  reflection  term. 


The  complete  evaluation  of  the  field  quantities  begins  with  the  twofold  integral  form  of  the 
dyadic  Green’s  function.  A  change  of  variable  is  applied  to  obtain  the  integration  in  terms  of  kp 
and  Bessel  functions.  In  integral  form,  the  resulting  expressions — for  which  exact  closed-form 
solutions  do  not  exist — of  the  x,  y,  and  z  components  of  the  scattered  field  are  given  in 
appendix  B.  Since  numerical  computation  of  these  expressions  is  slow  and  formidable  due  to  the 
presence  of  singularities  and  the  highly  oscillatory  behavior  of  the  Bessel  function  at  large 
distance  p,  approximate  asymptotic  solutions  are  sought.  In  a  standard  procedure,  Hankel’s 
functions  of  the  first  and  second  kind  are  used  in  extending  the  limits  of  integration  to  negative 
and  positive  infinity  (in  /cp-plane).  By  transforming  the  integration  over  to  the  complex  w-plane 
by  the  change  of  variable  kp  =  k0  sin  w,  and  then  substituting  the  asymptotic  form  of  the  Hanker  s 
function  for  large  arguments,  the  integrand  can  be  written  in  the  following  form — to  which  the 
standard  method  of  steepest  descents  can  be  applied: 


E-j  (p,  <f>,  z) 


zKJhLk 

1  6tt 
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where  subscripts  i,  j  =  x,  y,  or  z  and  q( w)  =  ikaR  cos{w-ws).  The  saddle  point  at  w  =  ws  is  defined 
by  q  ’(Wj)  =  0;  therefore 
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sinw5  =  —, 
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cos  w.  =  - 


z  +  z 
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R2=p2+(z  +  z')2. 
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After  some  manipulations  and  following  the  procedure  of  saddle  point  integration  provided  by 
Felsen  and  Marcuvitz  (48),  the  j  component  of  the  saddle  point  contribution  up  to  second  order 
for  a  dipole  pointed  in  the  i  direction  can  be  written  as 
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where  the  two  terms  represent  the  first  and  second  saddle  point  expansion  terms.  The  first  saddle 
point  expansion  term  is  the  geometrical-optics  reflection  term  and  the  second  expansion  term  can 
be  thought  as  a  correction  term  that  becomes  dominant  when  the  source  and  observation  points 
are  close  to  the  interface.  This  correction  term  is  also  commonly  known  as  the  Norton  wave  (49). 
After  some  lengthy  algebraic  manipulations,  expressions  for  fij(w)  of  (A-9)  can  be  shown  to  be 
the  following: 


2  - 

fzx= - sin 2  wcos  w  Rjm  (w) cos  (j>, 

k„ 


f  =  — —  sin 2  wcos  w Rj'm  ( w) sin <j), 
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k„ 


fzz  =  7—  sin 2  w  Rjm  (w), 


(A- 10) 
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fxz=  —  sin 2  wcos  w  Rj1m  (w)  cos  (j), 
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The  mixed  reflection  coefficient,  which  is  required  for  an  v  or  y-directed  dipole,  is  defined  as 


R+(w)  ^iW±^w<(w)_ 

ka  cos  w  ka 

The  TE  and  TM  reflection  coefficients  are  of  the  standard  form: 
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(A- 18) 

(A- 19) 


The  appearance  of  branch  points  at  wb  =  ±sin  ‘  entails  the  use  of  a  two-sheeted  Riemann 

surface  in  representing  the  w-plane;  in  order  to  satisfy  the  radiation  condition,  the  path  of 
integration  of  (A-5)  is  restricted  to  lie  on  the  upper  sheet  ( im  -  sin  2w  >o).  The  branch  cut 

defined  by  wb  =  sin-1  may  be  crossed  by  the  steepest  descent  path  and  its  contribution  to  the 
integral  can  be  included  by  adding  the  following  to  the  saddle  point  contribution  (48): 
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It  can  be  readily  verified  (48,  50)  that  the  poles  of  the  integrand  in  equation  A-5  are  not 
intercepted  as  the  original  integration  path  deforms  to  the  steepest  descent  path.  Depending  on 
the  permittivity  of  the  lower  medium  (cr|),  these  poles  (the  zeros  of  the  denominator  of  the 
reflection  coefficients  in  equations  A- 18  and  A- 19)  may  move  into  the  vicinity  of  the  saddle 
point;  thus,  the  poles  may  come  into  effect  indirectly  and  must  be  taken  into  account  by  carrying 
out  a  modified  saddle  point  integration  method.  Simulation  results,  however,  show  that  the 
ordinary  saddle  point  integration  as  used  here  is  sufficient  for  predicting  the  field  contributions  in 
the  far  field.  The  total  field  is  then  the  sum  of  the  direct,  saddle  point,  and  branch  cut 
contributions: 


Ev=Edv+E*+U(w,-0b)Ebs  (A-21) 

where  Ob  =  Re( w/,)  -  cos~' sech[\m(w h)\,  andU(')  is  the  Heaviside  step  function.  Calculation  of 
equation  A-20  can  be  quite  involved,  but  note  that  the  branch  cut  contribution  decays 
exponentially  with  p  since  the  wave  number  k\  is  complex;  therefore,  at  large  p,  the  branch  cut 
contribution  can  be  considered  as  negligible. 

For  a  two-layer  medium,  to  account  for  higher-order  reflections  that  are  transmitted  into  and  then 
emerging  from  the  dielectric  layer,  the  total  reflection  coefficient  is  modified  as  follows  (5i): 

Rp(w)  =  R°p  +  Xrf01(^2)'  ^yiTwei2udK^~5^ _  (A-22) 


Using  the  relations 


Rmn  =  _Rnm  ^  jmn  =  |  +  R,»n  ^ 


Equation  A-22  can  be  re-written  as 
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The  subscribe  “P”  is  a  place  holder  for  the  polarization — either  TE  or  TM.  The  saddle  point 
contribution  for  each  higher-order  reflection  (u  =  1,2,  . .  .qo)  in  the  summation  series  of 
equation  A-22  can  be  evaluated  by  using  equation  A-9  after  replacing  the  reflection  coefficient 
term  in  fjj(w)  with  t°  1  (R'P2 )“  (R'p  )u~lTp °ei2udk°  A«-si”2  w .  Note  that  the  higher-order  reflections  lead  to 
higher-order  poles  for  the  terms  in  the  summation  series  of  equation  A-22.  Although  they  are  still 
located  outside  of  the  steepest  descent  path,  these  poles  are  close  enough  that — at  higher 
orders — they  demand  a  more  refined  saddle  point  integration  method  even  when  p  is  relatively 
large.  In  order  to  circumvent  this  difficulty,  calculations  are  made  by  following  the  normal  mode 
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approach  in  which  the  total  aggregate  reflection  coefficient  in  equation  A-24  is  inserted  into 
equations  A-10-A-16  as  a  whole  rather  than  as  individuals  in  a  series.  This  approach  leads  to 
satisfactory  results  in  the  far  field  even  if  we  are  ignoring  the  pole  contributions — which  are  now 
supplied  by  the  zeros  of  the  denominator  in  equation  A-24  and  must  be  located  through 
numerical  methods  before  their  contributions  can  be  included  using  a  standard  technique  (48). 
These  pole  contributions,  depending  on  their  locations  on  the  complex  plane,  represent  either 
distinct  surface-wave  modes  or  leaky  modes — both  of  which  become  less  significant  as  the 
distance  between  source  and  observation  points  increases  (48,  52).  Also  note  that  in  the  two- 
layer  case  there  are  branch  points  only  at  wb  =  +sin_1  since  the  expression  in  equation  A-24 

can  be  shown  to  be  an  even  function  of  J£rl  -sin2  w  ■  The  branch  point  that  can  be  intercepted  is  at 
wh  =  sin-1  ,  but  since  medium  2  (the  ground  layer)  is  highly  lossy,  the  branch  cut  contribution 
falls  off  asymptotically  according  to  ~  eiklP  and,  hence,  rapidly  becomes  much  smaller  than  the 
algebraically  decaying  saddle  point  contribution. 

A-2.  Case  2  (-d  <  z’  <  0  and  z  >  0) 

Exact  formulations  for  this  case  have  been  derived  and  can  be  found  in  appendix  B.  Proceeding 
through  the  same  procedure  as  before,  the  method  of  steepest  descents  is  employed  in  obtaining 
the  integral  after  transforming  the  integration  to  the  w-plane.  Multiple  reflections  occurring 
within  the  dielectric  layer  can  be  accommodated  by  defining  the  total  transmission  coefficient  as 
the  following: 
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The  definition  of  equation  A-26  differs  from  that  of  equation  A-25  in  that  the  former  expression 
applies  to  waves  containing  an  initial  bounce  off  the  dielectric  layer/ground  interface.  It  can  be 
shown  that  the  transmitted  field  Ey  for  each  order  u  can  be  written  as 
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where 
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and  the  saddle  point  can  be  approximately  defined  by 
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The  set  of  fij,u(w)  for  an  arbitrary  oriented  dipole  can  be  shown  to  take  the  following  forms: 
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The  designations  “up”  and  “down”  differentiate  waves  that  are  initially  propagating  upward  and 
downward  from  the  source.  The  mixed  transmission  coefficients  have  been  written  in  the 
following  form: 
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In  the  formulations  above,  Tppu  and  r/7"  are  the  individual  terms  in  the  infinite  series  of  equations 

A-25  and  A- 26,  respectively.  Instead  of  computing  each  order  of  transmission  separately,  the 
saddle  point  evaluation  of  equation  A-27  can  be  carried  out  using  equation  A-9  as  before  by 
means  of  the  normal  mode  approach — in  which  now  the  aggregate  transmission  coefficients  in 
equations  A-25  and  A-26  are  substituted  into  fi/w)  as  one  term.  The  difficulty  mainly  lies  in 
taking  the  double  derivative  of  fijiw),  but  this  can  be  overcome  with  the  help  of  a  symbolic  math 
software.  When  the  source  and  observation  points  are  located  in  the  vicinity  of  the  interface,  it  is 
seen  that  both  the  first  and  second  term  in  the  saddle  point  expansion  fall  off  as  l/p2;  thus,  as  in 
Case  1,  both  expansion  terms  are  necessary  for  accurate  representation  of  the  total  field. 

It  is  seen  that  the  sign  of  the  term  j£rl  -sin2  w  has  no  effect  on  the  final  result  in  computing  the 

integral  in  equation  A-27;  therefore,  once  again,  the  function  fijiw)  is  an  even  function  of 
^srl  -  sin’  w  and  the  only  branch  points  on  the  complex  vv-plane  are  attributed  to  the  term 

Jer2  -  sin 2  w .  (As  a  matter  of  fact,  for  general  stratified  media  problems,  in  the  normal  mode 

approach,  the  branch  points  on  the  kp  plane  are  supplied  only  by  the  first  and  last  layer  (53).  The 
branch  points  of  the  first  layer  can  be  eliminated — as  it  has  been  done  here — by  translating  the 
calculation  onto  the  vv-plane  after  the  change  of  variable  kp  =  k0  sin  w.  Thus,  on  the  vv-plane,  the 
only  branch  points  remaining  are  due  to  the  last  layer  in  the  stratification.)  Similar  to  case  1, 
calculations  are  much  simplified  by  ignoring  contributions  from  the  branch  cut  and  the  poles 
(which  are  provided  by  the  zeros  of  the  denominator  of  equations  A-25  and  A-26).  In  the  far 
field,  this  claim  is  justified  since  the  saddle  point  contribution  becomes  the  only  dominant  field 
component. 

A-3.  Case  3  (-d  <  z’  <  0  and  -d  <  z  <  0) 

As  it  has  been  shown  in  reference  45,  when  both  the  source  and  observation  points  are  located 
inside  medium  2,  the  direct  and  saddle  point  contributions  decay  exponentially  as  a  function  of 
the  radial  distance  p  since  their  propagation  takes  place  in  a  lossy  medium.  The  dominant 
contribution  to  the  total  field,  as  it  turns  out,  comes  from  a  branch  cut  contribution.  Although  a 
modal  analysis  can  be  used,  a  ray  tracing  approach  provides  a  more  insightful  interpretation  for 
this  case.  In  the  ray  tracing  approach,  the  effective  reflection  coefficient  is  expanded  as  a  series 
before  asymptotic  evaluation  is  carried  out  separately  for  each  term.  In  such  an  approach,  both 
medium  1  and  medium  2  would  furnish  branch  cut  contributions  on  the  vv-plane  since  now — for 
each  individual  order  of  reflection — the  integrand  in  the  field  integral  is  no  longer  an  even 
function  of  J£rl  -  sin2  w .  Discarding  the  branch  cut  contribution  arising  from  the  branch  point  at 
wb  =  sin  1  yfsf  for  the  reason  mentioned  earlier,  the  only  relevant  and  significant  contribution,  in 
the  far  field,  comes  from  the  branch  cut  contribution  due  to  the  branch  point  at  wb  =  sin  1  fff  ■ 

This  branch  cut  contribution,  which  now  no  longer  undergoes  exponential  decay,  can  be 
interpreted  as  a  wave  from  the  source  that  radiates  upward  to  the  dielectric/air  interface  at  critical 
angle  and  then  propagates  along  the  interface  (in  air)  before  reaching  the  observation  point  at 
critical  angle  again.  Since  the  majority  of  the  propagation  takes  place  in  air,  this  field 
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component — which  has  been  labeled  as  the  “lateral”  wave — does  not  suffer  the  large  path  loss 
experienced  by  the  direct  and  saddle  point  contributions.  Through  standard  branch  cut 
integration  techniques,  the  three  lowest  orders  of  lateral  waves  have  been  derived  in  reference  45 
and  are  expressed  in  matrix  form  as  the  following  (z/  >  z): 
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where  a  and  are  symmetric  dyads  defined  in  reference  45.  The  first  term  in  equation  A-42 
represents  the  direct  lateral  wave  contribution;  the  second  term  is  the  contribution  of  the  lateral 
wave  generated  from  the  image  of  the  source  in  the  ground  plane;  and  the  third  term  is  the 
contribution  of  the  direct  lateral  wave  that  has  been  reflected  from  the  ground  plane  before 
reaching  the  observation  point.  It  is  seen  that  the  lateral  waves,  and  hence  the  total  field, 
decreases  as  1/p2 — which  is  the  same  asymptotic  behavior  observed  for  the  first  two  cases.  For 
further  details  on  the  derivation  and  verification  of  equation  A-42,  the  reader  is  referred  to 
reference  45. 


Although  the  transmitter  has  been  restricted  to  an  electric  dipole  throughout  this  work,  extension 
to  an  arbitrary  radiating  source  can  be  made  by  noting  that  the  asymptotic  form  of  the  Green’s 
function  for  each  of  the  three  cases  discussed  is  related  to  the  derived  expressions  for  the  field 
components  by  a  simple  constant: 


G(r',r ) 
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Therefore,  once  the  current  distribution  of  an  arbitrary  source  is  known,  an  approximation  to  the 
far-field  pattern  can  be  easily  computed.  Also,  upon  application  of  the  reciprocity  principle,  it  is 
straightforward  to  verify  that 

G{7,7')  =  G(r',7f  (A-44) 

where  the  superscript  “T”  indicates  the  transpose  operation.  Simply  by  using  the  relation  stated 
in  equation  A-44,  the  formulations  can  be  extended  to  other  standard  configurations  not 
explicitly  treated  above — for  example,  the  dipole  is  located  in  air  while  the  observation  point  is 
located  inside  the  dielectric  layer. 
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A  comparison  of  the  asymptotic  solution  with  the  exact  solution  is  shown  in  figures  A-l  and 
A-2 — for  a  half-space  radiation  problem. 


Figure  A-l.  (a)  Magnitude  and  (b)  phase  of  the  electric  fields  of  a  horizontal 
electric  dipole  as  a  function  of  range.  The  dipole  is  located  on  the 
z-axis  at  a  2-m  height  above  the  ground  and  the  observation  point  is 
at  a  1.7-m  height  along  a  radial  line  at  §  =  30°.  Ground  properties: 
£r’  =  4,  crd  =  10  mS/m,  and/=  300  MHz. 
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Figure  A-2.  (a)  Magnitude  and  (b)  phase  of  the  electric  fields  of  a 

horizontal  electric  dipole  as  a  function  of  range.  The  dipole 
is  buried  in  the  ground  at  a  10-cm  depth  and  the  observation 
point  is  at  a  1.7-m  height  along  a  radial  line  at  §  =  30°.  Ground 
properties:  £r’  =  4,  erd  =  10  mS/m,  and/=  300  MHz. 
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Appendix  B.  Sommerfeld  Integrals  for  Electric  Current  Sources 


Exact  integral  representations  of  the  radiation  from  an  electric  dipole  located  in  the  presence  of  a 
two-layer  infinite  medium  are  presented  below.  The  details  of  the  derivation  can  be  found  in  the 
work  by  Liao  (24). 


B-l.  Integral  Solutions  for  Case  1  (z’  >  0  and  z  >  0) 

Exact  expressions  for  the  electric  field  of  a  dipole  radiating  in  the  presence  of  a  two-layer 
medium  are  derived  from  the  dyadic  Green’s  function.  If  the  arbitrarily  oriented  dipole  is  located 
on  the  z-axis  at  (0,  0,  z’)  with  the  current  moment  vector  i  j  =  io (/  i  +  /  y  +  l_z) ,  it  can  be  shown 

that — for  Case  1 — the  general  equation  for  the  scattered  electric  field  at  observation  point  (p,  </>, 
z)  with  z’  >  z  can  be  written  as 


E-j=~^r 


J 


(B-l) 


where  the  function  in  the  integrand  is  dependent  upon  the  dipole  orientation  and  the  field 
component  of  interest  (x,  y,  or  z): 


fzx (kp )  =  r  rtm  (kp )j  1  (kpp) cos</>, 


(B-2) 


fzy(kp)  =  '-^Rm(kp)Jl(kpP)sm^ 


(B-3) 


2k 

fJkp)  =  1-pR™(kp)J0(kpP), 


(B-4) 


fa.M  =  URJkJ.JJkp)  ±  R+  (k  )J  2{k  p)  cos  (2^)], 


(B-5) 


UyJkJ  =  kpR+(k)J2(kp)  sin(20. 


(B-6) 


f xz  (kp )  =  4^  Rm  (kp )./,  (kpp)cos</>, 
k„ 


(B-7) 


f yz 0 kp )  =  —p~ rtm  (kp)J i (kp p)sin</>. 
k. 


(B-8) 


The  total  reflection  coefficient  is  given  by 


oo  M  _ f  Z?01\2 -i  e>12 z'2klz4  p01  ,  jpl2  i2klzd 

u  _  dOI  ,  tOItIOV /d10\w-1/ r>i2j2klzd ^  _  Doi  ,  L1  VKp  )  J Kpe  '  _  np  +  Kp  e 

Rp=RP+TpTpyRp)  (Rpe  ^=Rp+—^ = 


(B-9) 


and  the  mixed  reflection  coefficient  is 
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R±  (kp )  =  —j~—  — 


RrE(kp)  RTM(kp)k0Z 


(B-10) 


where  the  simple  reflection  coefficients  for  a  wave  going  from  layer  m  to  layer  n  are  the 
following: 


If  _  If 

D  ^ mz  nz 

^ TE  7  ,  7 

A:  +  k 


(B-ll) 


r,™.  _  ~ 

^TM  1  ,  ? 

£  A.  +  £  £ 


(B-12) 


The  wavenumber  is  represented  by  ku(u  =  0,  1,  or  2)  and 


k,„  =  Jk„2  -  k  2 . 


(B-13) 


The  direct  field  must  be  added  to  the  scattered  field  to  obtain  the  total  field.  If  z’  <  z,  the  relation 
in  equation  A-44  should  be  applied. 

B-2.  Integral  Solutions  for  Case  2  (-d  <  z’  <  0  and  z  >  0) 

When  the  dipole  is  inside  the  dielectric  layer  and  the  observation  point  is  in  air,  the  derivation  for 
the  total  field  is  more  complicated  but  is  not  much  different  from  that  of  the  first  case.  Simple 
manipulations  of  the  dyadic  Green’s  function  give  the  transmitted  field  as 


COjuJJ , 


\f,j(kp)elk°’zdkp 


(B-14) 


where  the  set  of  Mkp)  functions  can  be  shown  to  be 


-  2  ik  k 

/„(*,)  =  -  i  °-[TZe,k'^  +  T^Vk'^ ]7, (kpp) cos^, 

k„  k. 


(B-15) 


fzy(kp)  =  - 


■2  ik„  k 


£ — -\Tj^eik'^-  +  ]J,  (kp)  sin^ 


(B-16) 


/ zz  (kp )  =  — PT  ra  +  TlTe*"  {2d+°  Vo  (kpPX 

k,k„ 


(B-17) 


/„,*(*„)  =  kp{\T?e**tt'  +Tdowne,k^a+z)]J0(kpp) 

±[TuPe*l»  +T*owneiku&d+z')]J2(kpp)cos(2</>)}, 


(B-18) 


/*.„  (kp )  =  kp  {r"V^ 1  +  Trne'k^a+Z  >  }J2  (kpP) sin(2^). 


(B-19) 


fxz  (kp )  =  — ]7,  (kpp)  cos^, 


(B-20) 


fyz  (kp )  =  — f  —  [“Tj-'SJe*1'12'1  +  (£  p)sin^. 


(B-21) 
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The  various  transmission  coefficients  are  defined  below: 


TZ(k  )  TZ(k  )k 

rrupfj  \  _  te  v'V'  _u  tm  ^p'^oz 

*(')=  ^  K 


Tfwn(kp)  =  - 


"(K)  KTXk^ 


Tpp =Tp°YJ(K°Rp2ei2kudy =- 

11=0  1 


01  pi  2  J2kud  ’ 


+  RyRl;e 


Td™n  =  =  - 

u—0  1 


7)12/7-110 

iVp  ip 


+  RfR?eal'’f‘ 


(B-22) 

(B-23) 

(B-24) 

(B-25) 
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Appendix  C.  Traveling- Wave  Amplifier  Design 


Detailed  procedure  for  the  design  of  the  metal  semiconductor  field  effect  transistor  (MESFET) 
distributed  amplifier  employed  in  the  simulation  example  of  section  3.1  is  outlined.  First,  a 
simplified  small-signal  model  is  derived,  then  preliminary  designs  are  achieved  by  following  the 
guidelines  reported  by  Beyer  et  al.  (54).  From  these  preliminary  designs,  it  is  shown  that  the 
optimum  amplifier  should  be  a  two-stage  network  (n  =  2). 


The  amplifier  is  based  on  the  TGF-4350EPU  transistor  from  TriQuint.  In  a  standard  procedure, 
the  small-signal  characteristics  of  the  device  must  be  obtained  at  the  beginning  of  the  design 
process.  The  5-parameters  of  the  transistor  are  calculated  using  the  full  small-signal  model 
(figure  C-l)  provided  by  TriQuint. 


Figure  C-l.  TriQuint  small-signal  model  for  transistor. 


Next,  a  simplified  small-signal  model  for  the  transistor  is  derived  using  a  circuit  simulator  and  is 
shown  in  figure  C-2. 
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Figure  C-2.  Simulation  with  simplified  small-signal  model. 


The  5-parameters  of  the  simplified  model  are  compared  to  those  of  the  original  complete  model 
in  figure  C-3.  It  can  be  seen  that  there  is  a  very  good  match  for  5n  and  52 1;  there  is  a  poor  match 
for  S22  since  the  original  model  contains  parasitics  that  cannot  be  included  in  the  simplified 
model;  also  note  that  the  simplified  model  is  unilateral. 
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Figure  C-3.  Comparison  between  full  small-signal  model  (blue  lines)  and  simplified  small-signal  model  (red 
lines). 


Using  the  values  of  the  components  in  the  simplified  model,  the  design  procedure  outlined  by 
Beyer  et  al.  is  followed.  The  initial  design  curves  are  plotted  in  figure  C-4.  For  a  50  Q  design,  the 
cutoff  frequency  can  be  calculated  with  the  following  relations  (54): 


*01  — 


tfcC. 


~  *02 


I 


rfc(Cds  +  Cp) 


=  50Q. 


(C-l) 


Since  Cgs  is  fixed,  fc  must  be  equal  to  10.55  GHz.  It  is  easily  seen  that,  to  equalize  the  impedance 
of  the  gate  and  drain  lines,  an  additional  capacitance,  Cp  =  0.5784  pF,  must  be  added  to  the  drain. 
The  required  inductances  on  the  gate  and  drain  lines  are  then  calculated  as  (54) 

fc=  /  =  ,  1  =  Ld  =1.5085  nH.  (C-2) 

Lg  C  gs  Ld  (Cds  +  Cp ) 

As  it  is  done  in  reference  54,  the  image  terminations  are  realized  by  using  m-derived  filter  half 
sections  (55)  at  both  ends  of  the  lines.  The  form  of  the  termination  is  shown  in  figure  C-5.  The 
values  of  the  inductors  and  capacitor  are 

Lm  =^  =  0.4526  nH;  (C-3) 

L(\-m 2) 

Lm2  =  - -  =  0.8045  nH  ;  (C-4) 

2m 
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mC„, 

Cm= - —  =  0.1810/?F  ; 


where,  as  shown  in  reference  55,  the  optimum  value  for  m  is  0.6. 

j-/c  =  23.6  GHz 


(C-5) 


fc  =  13  GHz 


K  =  0.40 
AT  =0.35 


Figure  C-4.  Design  curves. 
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Figure  C-5.  Line  termination. 

Once  the  small-signal  model  of  the  device  and  the  values  of  the  different  components  (Lg,  Ld,  Cp, 
Lm,  Lm2,  Cm)  in  the  amplifier  network  are  defined,  the  preliminary  design  is  carried  out  in  the 
circuit  simulator.  An  example  of  the  preliminary  design  using  lumped-elements  is  shown  in 
figure  C-6  (for  n  =  2).  The  schematics  of  the  distributed  designs — for  n  =  1,  2,  3,  4 — are  also 
derived  (figure  C-7  shows  the  design  for  n  =  2).  In  each  of  these  preliminary  distributed  designs, 
inductors  are  replaced  by  short  transmission  lines  with  parameters  Za=  100  Q  and  L  =  ZJAgf; 
and  lumped-element  capacitors  are  still  used  (but  can  be  replaced  with  MIM  capacitors  later  in 
the  design  process).  Bias  is  done  with  an  ideal  DC-feed  inductor  and  DC-block  capacitors  are 
added.  Also,  note  that  the  actual  (not  simplified)  small-signal  model  is  used  for  the  transistor  in 
these  designs — all  of  which  are  over  a  gallium  arsenide  (GaAs)  substrate.  After  making  the 
above  changes,  the  component  values  of  each  circuit  are  tuned  to  maximize  the  gain  while 
keeping  its  variation  less  than  1  dB  over  the  entire  bandwidth.  It  is  seen  that  a  bandwidth  of 
10  GHz  can  be  realized  with  a  design  using  n  =  2.  With/C  fixed,  a  larger  number  of  devices 
increases  the  gain  but  reduces  the  bandwidth,  while  a  smaller  number  of  devices  leads  to  a 
smaller  gain  but  wider  bandwidth — this  is  consistent  with  the  trends  displayed  in  figure  C-4. 
Consequently,  from  a  survey  of  these  preliminary  designs,  with  the  constraint  that  the  bandwidth 
equals  10  GHz,  it  is  concluded  that  the  optimum  number  of  devices  is  n  =  2. 
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Figure  C-6.  (a)  An  example  layout  for  n  =  2  using  lumped-elements,  (b)  component  values,  and  (c)  simulation 
results  of  lumped-element-based  design. 
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Figure  C-7.  Schematic  and  gain  simulation  for  preliminary  design  with  n  =  2. 
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List  of  Symbols,  Abbreviations,  and  Acronyms 


3-D 

three-dimensional 

FDTD 

finite-difference  time-domain 

MESFET 

metal  semiconductor  field  effect  transistor 

MoM 

method-of-moments 

MPIE 

mixed  potential  integral  equation 

PCBs 

printed  circuit  boards 

RF 

radio  frequency 

SPICE 

Simulation  Program  with  Integrated  Circuit  Emphasis 

TE 

transverse  electric 

TM 

transverse  magnetic 
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